home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / yamp16.zip / VIRTDISK.CPP < prev    next >
C/C++ Source or Header  |  1993-01-11  |  28KB  |  1,197 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. // #include "virt.h"
  21. #include "alloc.h"
  22. //////////////////////////////////////////// string functions
  23.  
  24. strtype::strtype(strtype &str)     // copy constructor
  25. { int len = MAXCHARS;
  26.     len = strlen(str.name);
  27.     len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
  28.     strncpy(name, str.name, len);
  29.     name[len] = '\0';
  30. }
  31. strtype::strtype(char *str)
  32. { int len = MAXCHARS;
  33.     len = strlen(str);
  34.     len = (len >= MAXCHARS) ? MAXCHARS - 2 : len;
  35.     strncpy(name, str, len);
  36.     name[len] = '\0';
  37. }
  38. strtype::strtype(void)
  39. {
  40.     name[0] = '\0';
  41. }
  42.  
  43. strtype strtype::operator + (strtype str)
  44. {
  45.  
  46.     int len1, len2;
  47.     len1 = strlen(name);
  48.     len2 = strlen(str.name);
  49.     int len =(len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
  50.     len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
  51.     strncat(name, str.name, len);
  52.     name[len + len1] = '\0';
  53.     return *this;
  54. }
  55. strtype strtype::operator + (const char *str)
  56. {
  57.     int len1 = strlen(name);
  58.     int len2 = strlen(str);
  59.     int len = (len2 + len1 >= MAXCHARS) ? MAXCHARS - len1 - 3 : len2;
  60.     len = (len < 0 || len >= MAXCHARS - 1) ? 0 : len;
  61.     len = (len >= MAXCHARS) ? 0 : len;
  62.     strncat(name, str, len);
  63.     name[len + len1] = '\0';
  64.     return *this;
  65. }
  66. strtype strtype::operator = (strtype str)
  67. {
  68.     int len = (strlen(str.name) >= MAXCHARS) ? MAXCHARS - 3
  69.             : strlen(str.name);
  70.     strncpy(name, str.name, len);
  71.     name[len] = '\0';
  72.     return *this;
  73. }
  74. strtype strtype::operator = (const char *str)
  75. {
  76.     int len = (strlen(str) >= MAXCHARS) ? MAXCHARS - 2 : strlen(str);
  77.     strncpy(name, str, len);
  78.     name[len] = '\0';
  79.     return *this;
  80. }
  81.  
  82. //////////////////////////////////// virtual memory allocation functions
  83. // Virtual memory routines used by permission
  84. // Portions of this code are (c) 1991 by Allen I. Holub and are used by
  85. // permission
  86. // The code is found in Programmer's Journal volume 8.5. (1990)
  87. //
  88.  
  89.  
  90. //#define DEBUG 1
  91.  
  92. #ifdef DEBUG
  93. #define D(x) x
  94. #define ND(x)
  95. #else
  96. #define D(x)
  97. #define ND(x) x
  98. #endif
  99.  
  100. #define READ 1
  101. #define WRITE 0
  102. #define TNAME "$$$$vmem.tmp"
  103.  
  104. #define BLK_SIZE 512
  105. #define SIGNATURE 0x5a5a
  106.  
  107. static char Vname[128];
  108. static unsigned Filesize = 0;
  109. static int Vfd = -1;
  110. static int nobj = 0;
  111.  
  112. static hdr *Freelist = NULL;
  113.  
  114. static hdr *new_hdr(unsigned nblocks, unsigned offset);
  115. static unsigned enlarge(unsigned blocks_reqd);
  116. static hdr *add_vmem(long num_ele, int ele_size, hdr *header);
  117. static int load_block(hdr *p, unsigned req_block);
  118. static int block_io(int doread, unsigned block, char *buf);
  119.  
  120. /*--------------------------------------------*/
  121. void pheader(hdr *p)
  122. {
  123.     char junk[12] = " <-Freelist";
  124.     char s[13] = " ";
  125.  
  126.     if (!p) printf("NULL\n");
  127.     else {
  128.         if (p == Freelist) strcpy(s, junk);
  129.         printf("hdr  0x%lX :active %1d offset = %u, nblocks = %u, next=0x%lX %s\n",
  130.             p, p->active, p-> offset, p-> nblocks, p-> next.h, s);
  131.     }
  132. }
  133. void pfreelist(void)
  134. {
  135.     hdr *p = NULL;
  136.     if (!(p = Freelist)) printf("Freelist is empty.\n");
  137.     else {
  138.         printf("Freelist is:\n");
  139.         do {
  140.             pheader(p);
  141.         } while ((p = p->next.h) != Freelist);
  142.     }
  143. }
  144.  
  145. /*-------------------------------------*/
  146. int vopen(void)
  147. {
  148.     char *env;
  149.     int len;
  150.     if (!(env = getenv("TMP")))
  151.         strcpy(Vname, TNAME);
  152.     else {
  153.         if ((len = strlen(env)) > sizeof (Vname) - 16) {
  154.             errno = E2BIG;
  155.             return 0;
  156.         }
  157.         sprintf(Vname,
  158.                 (len >= 2 && env[len - 2] != '/' && env[len - 2] !=    '\\')?
  159.             "%s/%s" : "%s%s", env, TNAME);
  160.     }
  161.     if ((Vfd = creat(Vname, S_IREAD | S_IWRITE)) == -1) return 0;
  162.     close(Vfd);
  163.     if ((Vfd = open(Vname, O_RDWR | O_BINARY)) == -1) return 0;
  164.     return 1;
  165. }
  166. /*---------------------------------*/
  167. int vclose(void)
  168. {
  169.     hdr *p = NULL, *newp = NULL;
  170.     if (Vfd == -1) {
  171.         errno = EBADF;
  172.         return 0;
  173.     }
  174.     else {
  175.         if (close(Vfd) == -1) return 0;
  176.         Vfd = -1;
  177.         ND(if (unlink(Vname) == -1) return 0;)
  178.         if (Freelist) {
  179.             p = Freelist;
  180.             do {
  181.                 newp = p->next.h;
  182.                 free(p);
  183.             } while ((p = newp) != Freelist);
  184.             Freelist = NULL;
  185.         }
  186.     }
  187.     return 1;
  188. }
  189. /*----------------------------------*/
  190. hdr *vmalloc(long num_ele, size_t ele_size)
  191. {
  192.     hdr *prev, *cur, *newhdr, *eof_block;
  193.     unsigned ele_per_block;
  194.     unsigned blocks_reqd;
  195.  
  196.     D(printf("Entering vmalloc(%ld,%u). ", num_ele, ele_size);)
  197.     D(pfreelist();)
  198.     ele_per_block = BLK_SIZE / ele_size;
  199.     blocks_reqd = num_ele / ele_per_block + (num_ele % ele_per_block != 0);
  200.  
  201.     if (Freelist) {
  202.         prev = Freelist;
  203.         cur = Freelist->next.h;
  204.         eof_block = NULL;
  205.  
  206.         while (1) {
  207.             do {
  208.                 if (cur-> offset + cur-> nblocks == Filesize) eof_block = cur;
  209.                 if (cur-> nblocks >= blocks_reqd) {
  210.                     Freelist = prev;
  211.                     if (cur-> nblocks > blocks_reqd) {
  212.                         cur->nblocks -= blocks_reqd;
  213.                         if (newhdr = new_hdr(blocks_reqd, cur->offset + cur->nblocks))
  214.                             newhdr = add_vmem(num_ele, ele_size, newhdr);
  215.                         goto end;
  216.                     }
  217.                     else {
  218.                         if (cur-> next.h == cur) Freelist = NULL;
  219.                         else prev->next.h = cur->next.h;
  220.                         newhdr = add_vmem(num_ele, ele_size, cur);
  221.                         goto end;
  222.                     }
  223.                 }
  224.                 prev = cur;
  225.                 cur = cur->next.h;
  226.             } while (prev != Freelist);
  227.  
  228.             if (!eof_block) break;
  229.             else if (!enlarge(blocks_reqd - eof_block->nblocks)) {
  230.                 newhdr = NULL;
  231.                 goto end;
  232.             }
  233.             else eof_block->nblocks = blocks_reqd;
  234.         }     /* end while(1) */
  235.     }     /* end  if(freelist) */
  236.  
  237.     if (!(newhdr = new_hdr(blocks_reqd, Filesize))) goto end;
  238.     if (!enlarge(blocks_reqd)) {
  239.         free(newhdr);
  240.         newhdr = NULL;
  241.     }
  242.     newhdr = add_vmem(num_ele, ele_size, newhdr);
  243.     end :
  244.     if (newhdr != NULL) {
  245.         newhdr->active = 1;
  246.         newhdr->nrefs = 1;
  247.     }
  248.     D(printf("leaving vmalloc(%ld,%u), Return:\n", num_ele, ele_size);)
  249.     D(pheader(newhdr);)
  250.     D(pfreelist();)
  251.     D(printf("-----------------------------\n");)
  252.     return newhdr;
  253. }
  254. /*---------------------------------------*/
  255. hdr *new_hdr(unsigned nblocks, unsigned offset)
  256. {
  257.     hdr *newhdr = NULL;
  258.     if (!(newhdr = (hdr *) malloc(sizeof (hdr)))) {
  259.         errno = ENOMEM;
  260.         return NULL;
  261.     }
  262.     newhdr->nblocks = nblocks;
  263.     newhdr->offset = offset;
  264.     return newhdr;
  265. }
  266. /*---------------------------------------*/
  267. unsigned enlarge(unsigned blocks_reqd)
  268. {
  269.     unsigned size_in_blocks;
  270.     long real_size, position;
  271.  
  272.     position = ((long) blocks_reqd) * ((long) BLK_SIZE) - 1L;
  273.  
  274.     if ((real_size = lseek(Vfd, position, SEEK_END)) == -1) {
  275.         errno = ENOMEM;
  276.         return 0L;
  277.     }
  278.     if (write(Vfd, "", 1) == -1)
  279.         return 0L;
  280.  
  281.     ++ real_size;
  282.     return ((real_size / BLK_SIZE > (unsigned) ~0) ? 0 :
  283.         (Filesize = (unsigned) (real_size / BLK_SIZE)));
  284. }
  285. /*---------------------------------------*/
  286. hdr *add_vmem(long num_ele, int ele_size, hdr *header)
  287. {
  288.     vmem *p = NULL;
  289.  
  290.     if (!(p = (vmem *) malloc(sizeof (vmem)))) {
  291.         errno = ENOMEM;
  292.         return NULL;
  293.     }
  294.     else if (!(p->buf = (char *) malloc(BLK_SIZE))) {
  295.         errno = ENOMEM;
  296.         return NULL;
  297.     }
  298.     else {
  299.         header->next.v = p;
  300.         p->signature = SIGNATURE;
  301.         p->dirty = 0;
  302.         p->ele_size = ele_size;
  303.         p->ele_per_blk = BLK_SIZE / ele_size;
  304.         p->num_ele = num_ele;
  305.         p->cblock = 0;
  306.         return header;
  307.     }
  308. }
  309. /*-------------------------------------*/
  310. int vfree(hdr *p)
  311. {
  312.     int rval = 1;
  313.     hdr *cur, *h;
  314.     vmem *v;
  315.  
  316.     v = p->next.v;
  317.     D(printf("Entering vfree(0x%lX).\n", p);)
  318.     D(printf("disposing: "); pheader(p);)
  319.     D(pfreelist();)
  320.  
  321.     if (Vfd == -1) {     /* VMS system is not open */
  322.         D(printf("VFREE: VMS system closed\n");)
  323.         goto end;
  324.     }
  325.  
  326.     p->nrefs = 0;
  327.     if (v-> signature == SIGNATURE) v-> signature = 0;
  328.     else {
  329.         errno = EBADF;
  330.         rval = 0;
  331.         goto end;
  332.     }
  333.  
  334.     p->active = 0;
  335.     free(v->buf);
  336.     free(v);
  337.     if (!Freelist) {
  338.         Freelist = p->next.h = p;
  339.         goto end;
  340.     }
  341.  
  342.     for (cur = Freelist; 1; cur = cur-> next.h) {
  343.         if ((cur->offset < p->offset && p->offset < cur->next.h->offset)
  344.             || (cur-> offset >= cur-> next.h->offset &&
  345.                 (cur-> offset < p-> offset ||
  346.                     p->offset < cur-> next.h->offset)))
  347.             break;
  348.     }
  349.     if (cur == cur-> next.h) {
  350.         if (p-> offset + p-> nblocks == cur-> offset) {
  351.             cur->offset = p-> offset;
  352.             cur->nblocks += p-> nblocks;
  353.             if (p == Freelist) Freelist = cur;
  354.             free(p);
  355.             goto end;
  356.         }
  357.         else if (cur-> offset + cur-> nblocks == p-> offset) {
  358.             cur->nblocks += p-> nblocks;
  359.             if (p == Freelist) Freelist = cur;
  360.             free(p);
  361.             goto end;
  362.         }
  363.     }
  364.     if (p-> offset + p-> nblocks == cur-> next.h->offset) {
  365.         p->nblocks += cur-> next.h->nblocks;
  366.         p->next.h = cur->next.h->next.h;
  367.         if (cur-> next.h == Freelist) Freelist = cur;
  368.         free(cur->next.h);
  369.     }
  370.     else
  371.         p->next.h = cur->next.h;
  372.  
  373.     if (cur-> offset + cur-> nblocks == p-> offset) {
  374.         cur->nblocks += p-> nblocks;
  375.         cur->next.h = p->next.h;
  376.         if (p == Freelist) Freelist = cur;
  377.         free(p);
  378.     }
  379.     else cur->next.h = p;
  380.  
  381.     end :
  382.     D(printf("leaving vfree( 0x%lX). Returning %d. ", p, rval);)
  383.     D(pfreelist();)
  384.     D(printf("-----------------------------\n");)
  385.     return rval;
  386. }
  387. /*-----------------------------------------------*/
  388. int block_io( int doread, unsigned block_num, char *buf)
  389.   {
  390.     size_t got;
  391.     long here;
  392.     here = lseek(Vfd, ((long)block_num * (long)BLK_SIZE) , SEEK_SET);
  393.     if (here == -1 ) return 0;
  394.     got = doread ? read (Vfd, buf, (size_t) BLK_SIZE)
  395.          : write(Vfd, buf, (size_t) BLK_SIZE);
  396.     if ( got == -1 ) return 0;
  397.     return (doread ? 1 : (got == BLK_SIZE));
  398.   }/*----------------------------------------------*/
  399. int load_block(hdr *p, unsigned req_block)
  400. {
  401.     vmem *v;
  402.     v = p->next.v;
  403.     if (req_block != v-> cblock) {
  404.         if (v-> dirty && !block_io(WRITE, p->offset + v->cblock, v->buf))
  405.             return 0;
  406.         if (!block_io(READ, p->offset + req_block, v->buf)) return 0;
  407.         v->dirty = 0;
  408.         v->cblock = req_block;
  409.     }
  410.     return 1;
  411. }
  412. /*-----------------------------------------*/
  413. void *vread(hdr *p, long index)
  414. {
  415.     unsigned block_num;
  416.     vmem *v;
  417.     v = p->next.v;
  418.  
  419.     if (index < 0 || index > v-> num_ele) {
  420.         errno = ERANGE;
  421.         return 0;
  422.     }
  423.     if (v-> signature != SIGNATURE || Vfd == -1) {
  424.         errno = EBADF;
  425.         return 0;
  426.     }
  427.     block_num = index / v->ele_per_blk;
  428.     if (!load_block(p, block_num)) return NULL;
  429.  
  430.     index -= block_num * v->ele_per_blk;
  431.     index *= v->ele_size;
  432.     return (v-> buf + index);
  433. }
  434.  
  435. /*---------------------------------------------*/
  436. void *vwrite(hdr *p, long index, char *thiss)
  437. {
  438.     int i;
  439.     char *bp = NULL;
  440.     void *startp = NULL;
  441.     vmem *v;
  442.     v = p->next.v;
  443.     if (startp = vread(p, index)) {
  444.         bp = (char *) startp;
  445.         for (i = v-> ele_size;-- i >= 0;) *bp++ = *thiss++;
  446.         v->dirty = 1;
  447.     }
  448.     return startp;
  449. }
  450. /*--------------------------------------------------*/
  451. long vele(hdr *p) { return p->next.v->num_ele; }
  452. void vdirty(hdr *p) { p->next.v->dirty = 1; }
  453. /*--------------------------------------------------*/
  454.  
  455.  
  456.  
  457.  
  458. ///////////////////////////////////////////// virtual double array
  459.  
  460.  
  461. static double val(vdoub &v)
  462. {
  463.     if (v.cur_ele) return *v.cur_ele;
  464.     cerr << "VMS: no previous index on vdoub (val)\n";
  465.     return 0;
  466. }
  467.  
  468. vdoub::vdoub(long array_size)
  469. {
  470.     if (nobj++ == 0) {
  471.         if (!vopen()) {
  472.             perror("VMS: unable to open virtural memory file");
  473.             exit(1);
  474.         }
  475.     }
  476.     cur_ele = NULL;
  477.     cur_index = 0;
  478.     signature = SIGNATURE;
  479.  
  480.     if (!(hdr = vmalloc(array_size, sizeof (double)))) {
  481.         perror("VMS allocation error 1");
  482.         exit(1);
  483.     }
  484. }
  485.  
  486. vdoub::vdoub(void)
  487. {
  488.     if (nobj++ == 0) {
  489.         if (!vopen()) {
  490.             perror("VMS: unable to open virtural memory file");
  491.             exit(1);
  492.         }
  493.     }
  494.     cur_ele = NULL;
  495.  
  496.     cur_index = 0;
  497.     signature = SIGNATURE;
  498.  
  499.     if (!(hdr = vmalloc(1L, sizeof (double)))) {
  500.         Dispatch->PrintStack();
  501.         perror("VMS allocation error 2");
  502.         exit(1);
  503.     }
  504. }
  505.  
  506. vdoub::vdoub(vdoub &src)
  507. {
  508.     if (nobj++ == 0) {
  509.         if (!vopen()) {
  510.             perror("VMS: unable to open virtural memory file");
  511.             exit(1);
  512.         }
  513.     }
  514.     double *p;
  515.     long numele = vele(src.hdr);
  516.     signature = SIGNATURE;
  517.     cur_ele = NULL;
  518.     cur_index = 0;
  519.  
  520.     if (!(hdr = vmalloc(numele, sizeof (double)))) {
  521.         perror(" VMS copy allocation error");
  522.         exit(1);
  523.     }
  524.     while (-- numele >= 0) {
  525.         if (!(p = (double *) vread(src.hdr, numele))) {
  526.             perror("VMS copy-access( read ) error");
  527.             exit(1);
  528.         }
  529.         if (!vwrite(hdr, numele, (char *) p)) {
  530.             perror("VMS copy-access( write ) error");
  531.             exit(1);
  532.         }
  533.     }
  534. }
  535.  
  536. vdoub::~vdoub(void)
  537. {
  538.     if (0 ==-- (hdr->nrefs)) {
  539.         if (!vfree(hdr)) {
  540.             perror("VMS deallocation error ");
  541.             exit(1);
  542.         }
  543.     }
  544.     signature = 0;     // make sure to change signature in destructor
  545.     if (-- nobj == 0) vclose();
  546. }
  547.  
  548. double vdoub::v(long index)
  549. {
  550.     if (!(cur_ele = (double *) vread(hdr, index))) {
  551.         cerr << "hdr, index: " << hdr << " " << index;
  552.         perror("VMS selection error 1");
  553.     }
  554.  
  555.     cur_index = index;
  556.     return *cur_ele;
  557. }
  558.  
  559. vdoub &vdoub::operator[](long index)
  560. {
  561.     if (!(cur_ele = (double *) vread(hdr, index))) {
  562.         cerr << "hdr, index: " << hdr << " " << index << " ";
  563.         perror("VMS selection error 2");
  564.         pfreelist();
  565.         exit(1);
  566.     }
  567.     cur_index = index;
  568.     return *this;
  569. }
  570.  
  571. double vdoub::operator = (double d)
  572. {
  573.     if (Garbage()) {
  574.         perror("VMS assignment error, source is Garbage(=)");
  575.         exit(1);
  576.     }
  577.     *cur_ele = d;
  578.     vdirty(hdr);
  579.     return d;
  580. }
  581.  
  582. double vdoub::operator = (vdoub &v)
  583.     // copy vector v to vector t
  584.     {
  585.         if (v.Garbage()) {
  586.             perror("VMS copy error, source is Garbage");
  587.             exit(1);
  588.         }
  589.         if (!Garbage()) vfree(hdr);     // replace hdr if t is active
  590.  
  591.         double *p;
  592.         long numele = vele(v.hdr);
  593.         signature = SIGNATURE;
  594.         *cur_ele = *v.cur_ele;
  595.         cur_index = v.cur_index;
  596.  
  597.         if (!(hdr = vmalloc(numele, sizeof (double)))) {
  598.             perror(" VMS copy allocation error");
  599.             exit(1);
  600.         }
  601.         while (-- numele >= 0) {
  602.             if (!(p = (double *) vread(v.hdr, numele))) {
  603.                 perror("VMS copy-access( read ) error");
  604.                 exit(1);
  605.             }
  606.             if (!vwrite(hdr, numele, (char *) p)) {
  607.                 perror("VMS copy-access( write ) error");
  608.                 exit(1);
  609.             }
  610.         }
  611.         return *cur_ele;
  612. }
  613.  
  614.  
  615. //////////////////////////////////////////////////matrix stack functions
  616.  
  617. MStack::MStack(void)
  618. {
  619.     static int cnter = 0;
  620.     next = NULL;
  621.     stackloc = 0;
  622.     level = 0;
  623.     if (!cnter) {
  624.         VMatrix::Nameit("DISPATCH");
  625.         cnter = 1;
  626.         level = 1;
  627.     }
  628. }
  629.  
  630. void MStack::Inclevel(void)
  631. {
  632.     Dispatch->level++;
  633. }
  634.  
  635. void MStack::Declevel(void)
  636. {
  637.  
  638.     (Dispatch-> level)--;
  639.     if (Dispatch-> level < 0) {
  640.         printf("ERROR: LEVEL has been decremented too often\n");
  641.         exit(1);
  642.     }
  643. }
  644.  
  645. void VMatrix::NewReference(VMatrix &ROp)
  646. {
  647.  
  648.     cur_ele = NULL;
  649.     cur_index = 0;
  650.     signature = SIGNATURE;
  651.     r = ROp.r;
  652.     c = ROp.c;
  653.     // release the current header and share it with ROp.hdr
  654.     PurgeVectors();
  655.     hdr = ROp.hdr;
  656.     hdr->nrefs += 1;
  657. }
  658. void MStack::Push(VMatrix &ROp)
  659. {
  660.     ROp.Garbage("Push");
  661.     D(printf("Pushing: "); Dispatch->PrintStack();)
  662.     MStack *temp = new MStack;
  663.     temp->Nameit(ROp.Getname(ROp));
  664.     temp->NewReference(ROp);
  665.     temp->next = Dispatch-> next;
  666.     Dispatch->next = temp;
  667.     temp->level = Dispatch-> level;
  668.     temp->stackloc =++ (Dispatch->stackloc);
  669. }
  670.  
  671. VMatrix& MStack::Pop(void)
  672. {
  673.     D(printf("Popping: "); Dispatch->PrintStack();)
  674.  
  675.     if (Dispatch-> next && Dispatch-> stackloc) {
  676.         MStack *t = Dispatch->next;
  677.         Dispatch->next = Dispatch-> next-> next;
  678.         delete t;
  679.         (Dispatch-> stackloc)--;
  680.     }
  681.     else Nrerror("POP: Stack is empty, cant pop any more\n");
  682.     return *Dispatch;
  683. }
  684.  
  685. void MStack::Cleanstack(void)
  686. {
  687.     while (Dispatch-> next-> level >= Dispatch-> level
  688.         && Dispatch->next-> next)
  689.         Dispatch->Pop();
  690. }
  691. void MStack::PrintStack(void)
  692. {
  693.     MStack *temp = Dispatch;
  694.     int t = 1 + Dispatch->stackloc;
  695.  
  696.     while (t--) {
  697.         printf("\n Stack matrix :location level r c:%d %d %d %d: name: ",
  698.             t, temp->level, temp-> r, temp-> c);
  699.         temp->Showname();
  700.         temp = temp->next;
  701.     }
  702.     printf("\n\n");
  703. }
  704.  
  705. //////////////////////////////////////////////////////matrix class
  706.  
  707. VMatrix::VMatrix(void)
  708. {
  709.     r = c = 1;
  710.     curvecind = 0;
  711.     Nameit("t");
  712.     signature = SIGNATURE;
  713.     //    for( int i=1; i<=r; i++ ){
  714.     //      for(int j=1; j<=c; j++) M(i,j) = 0;
  715.     //    }
  716. }
  717.  
  718. VMatrix::VMatrix(const char *str, int rr, int cc) :
  719. vdoub(((long) ((long) rr) *((long) cc)))
  720. {
  721.     r = rr;
  722.     c = cc;
  723.     curvecind = 0;
  724.     Nameit(str);
  725.     signature = SIGNATURE;
  726.     //    for( int i=1; i<=r; i++ ){
  727.     //      for(int j=1; j<=c; j++) M(i,j) = 0;
  728.     //    }
  729. }
  730. VMatrix::VMatrix(VMatrix &ROp)     // copy constructor
  731. {
  732.     ROp.Garbage("Copy Constructor");
  733.     r = ROp.r;
  734.     c = ROp.c;
  735.     name = ROp.name;
  736.     curvecind = 0;
  737.     signature = SIGNATURE;
  738.     PurgeVectors();     // note a hdr is constructed in the
  739.     // vdoub constructor so it must be deleted
  740.     // first.
  741.     SetupVectors(r, c);
  742.     for (int i = 1; i <= r;++ i) {
  743.         for (int j = 1; j <= c;++ j) M(i, j) = ROp.m(i, j);
  744.     }
  745.     Dispatch->Cleanstack();
  746. }
  747. VMatrix::~VMatrix(void)
  748. {
  749.     signature = 0;
  750. }
  751.  
  752. //////////////////////////////////// matrix member functions
  753.  
  754. double VMatrix::Nrerror(const char *errormsg)
  755. {
  756.     double x;
  757.     printf("\nMATRIX ERROR: %s \nexiting to system\n", errormsg);
  758.     x = 2;
  759.     exit(1);
  760.     return x;
  761. }
  762. void VMatrix::Garbage(const char *errormsg) {
  763. #ifndef NO_CHECKING
  764.     int errorint = 0;
  765.     if (signature != SIGNATURE) errorint++;
  766.     if (errorint) {
  767.         Showname();
  768.         printf("\nFunction name: %s", errormsg);
  769.         Nrerror("Operating on Garbage");
  770.     }
  771. #endif
  772.     return;
  773. }
  774.  
  775. void VMatrix::SetupVectors(int rr, int cc)
  776. {
  777.     cur_ele = NULL;
  778.  
  779.     cur_index = 0;
  780.     signature = SIGNATURE;
  781.     r = rr;
  782.     c = cc;
  783.  
  784.     if (!(hdr = vmalloc(((long) ((long) r) *((long) c)), sizeof (double)))){
  785.         perror("VMS allocation error 3");
  786.         exit(1);
  787.     }
  788. }
  789.  
  790. void VMatrix::PurgeVectors(void)
  791. {
  792.     // frees the virtual vector
  793.     Garbage("PurgeVectors");
  794.     if (!vfree(hdr)) {
  795.         perror("VMS deallocation error ");
  796.         exit(1);
  797.     }
  798. }
  799. void VMatrix::DisplayMat(void)
  800. {
  801.     Garbage("DisplayMat");
  802.     int wid = Getwid();
  803.     int dec = Getdec();
  804.     printf("\n");
  805.     name.Showname();
  806.     printf("\n\n");
  807.     for (int i = 1; i <= r;++ i) {
  808.         for (int j = 1; j <= c;++ j) {
  809.             printf("%*.*lf ", wid, dec, m(i, j));
  810.         }
  811.         printf("\n");
  812.     }
  813.     printf("\n");
  814. }
  815. void VMatrix::InfoMat(void)
  816. {
  817.     Garbage("InfoMat");
  818.     printf("\n");
  819.     name.Showname();
  820.     printf("\nr c: %d %d\n", r, c);
  821. }
  822. void VMatrix::LoadMat(void)
  823. {
  824.     // interactive matrix loading;
  825.     char buffer[256] = { '\0' };
  826.  
  827.     Garbage("LoadMat");
  828.     printf(" Enter name: ");
  829.     cin >> buffer;
  830.     Nameit(buffer);
  831.  
  832.     double d;
  833.     for (int j = 1; j <= r;++ j) {
  834.         for (int k = 1; k <= c;++ k) {
  835.             cout << "Enter (" << j << "," << k << "): ";
  836.             cin >> d;
  837.             M(j, k) = d;
  838.         }
  839.     }
  840. }
  841.  
  842. void VMatrix::Replace(VMatrix &ROp)
  843. {
  844.     ROp.Garbage("Replace");
  845.     if (r != ROp.r || c != ROp.c) {
  846.         PurgeVectors();
  847.         SetupVectors(ROp.r, ROp.c);
  848.     }
  849.     name = ROp.name;
  850.     for (int i = 1; i <= r;++ i) {
  851.         for (int j = 1; j <= c;++ j) { M(i, j) = ROp.m(i, j);
  852.         }
  853.     }
  854. }
  855.  
  856. void VMatrix::operator = (VMatrix &ROp)
  857. {
  858.     Garbage("operator =");
  859.     ROp.Garbage("operator =");
  860.     D(printf("Equals "); Showname(); Dispatch->PrintStack();)
  861.     Replace(ROp);
  862.     Dispatch->Cleanstack();
  863. }
  864. ////////////////////// end matrix member functions
  865.  
  866. /////////////////// freind functions of matrix class
  867.  
  868. ///------------- addition
  869.  
  870. VMatrix& operator + (VMatrix &LOp, VMatrix &ROp)
  871. {
  872.     LOp.Garbage("operator +");
  873.     ROp.Garbage("operator +");
  874.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  875.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  876.         for (int i = 1; i <= LOp.r;++ i) {
  877.             for (int j = 1; j <= ROp.c;++ j) {
  878.                 temp->M(i, j) = LOp.m(i, j) + ROp.m(i, j);
  879.             }
  880.         }
  881.         temp->name = temp-> name + LOp.name + "+" + ROp.name + ")";
  882.         Dispatch->Push(*temp);
  883.         delete temp;
  884.     }
  885.     else ROp.Nrerror("Mismatched Matrix sizes in addition function\n");
  886.     return Dispatch->ReturnMat();
  887. }
  888. VMatrix& operator +(double scalar, VMatrix &ROp)
  889. {
  890.     ROp.Garbage("operator s+M");
  891.     strtype string;
  892.     char buffer[MAXCHARS] = { '\0' };
  893.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  894.         for (int i = 1; i <= ROp.r;++ i) {
  895.             for (int j = 1; j <= ROp.c;++ j) {
  896.                 temp->M(i, j) = scalar + ROp.m(i, j);
  897.             }
  898.         }
  899.         gcvt(scalar, NDECS, buffer);
  900.         string = buffer;
  901.         temp->name = temp-> name + string + "+" + ROp.name + ")";
  902.         Dispatch->Push(*temp);
  903.         delete temp;
  904.         return Dispatch->ReturnMat();
  905. }
  906. VMatrix& operator +(VMatrix &ROp, double scalar)
  907. {
  908.     ROp.Garbage("operator M+s");
  909.     strtype string;
  910.     char buffer[MAXCHARS] = { '\0' };
  911.  
  912.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  913.         for (int i = 1; i <= ROp.r;++ i) {
  914.             for (int j = 1; j <= ROp.c;++ j) {
  915.                 temp->M(i, j) = scalar + ROp.m(i, j);
  916.             }
  917.         }
  918.         gcvt(scalar, NDECS, buffer);
  919.         string = buffer;
  920.         temp->name = temp-> name + ROp.name + "+" + string + ")";
  921.         Dispatch->Push(*temp);
  922.         delete temp;
  923.         return Dispatch->ReturnMat();
  924. }
  925.  
  926. ////-------------subtraction
  927.  
  928. VMatrix& operator -(VMatrix &LOp, VMatrix &ROp)
  929. {
  930.     LOp.Garbage("operator M-M");
  931.     ROp.Garbage("operator M-M");
  932.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  933.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  934.         for (int i = 1; i <= LOp.r;++ i) {
  935.             for (int j = 1; j <= LOp.c;++ j) {
  936.                 temp->M(i, j) = LOp.m(i, j) - ROp.m(i, j);
  937.             }
  938.         }
  939.         temp->name = temp-> name + LOp.name + "-" + ROp.name + ")";
  940.         Dispatch->Push(*temp);
  941.         delete temp;
  942.     }
  943.     else Dispatch->Nrerror("Mismatched VMatrix sizes in subtraction function\n");
  944.     return Dispatch->ReturnMat();
  945. }
  946. VMatrix& operator -(double scalar, VMatrix &ROp)
  947. {
  948.     ROp.Garbage("operator s-M");
  949.     strtype string;
  950.     char buffer[MAXCHARS] = { '\0' };
  951.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  952.         for (int i = 1; i <= ROp.r;++ i) {
  953.             for (int j = 1; j <= ROp.c;++ j) {
  954.                 temp->M(i, j) = scalar - ROp.m(i, j);
  955.             }
  956.         }
  957.         gcvt(scalar, NDECS, buffer);
  958.         string = buffer;
  959.         temp->name = temp-> name + string + "-" + ROp.name + ")";
  960.         Dispatch->Push(*temp);
  961.         delete temp;
  962.         return Dispatch->ReturnMat();
  963. }
  964. VMatrix& operator -(VMatrix &ROp, double scalar)
  965. {
  966.     ROp.Garbage("operator M-s");
  967.     strtype string;
  968.     char buffer[MAXCHARS] = { '\0' };
  969.  
  970.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  971.         for (int i = 1; i <= ROp.r;++ i) {
  972.             for (int j = 1; j <= ROp.c;++ j) {
  973.                 temp->M(i, j) = ROp.m(i, j) - scalar;
  974.             }
  975.         }
  976.         gcvt(scalar, NDECS, buffer);
  977.         string = buffer;
  978.         temp->name = temp-> name + ROp.name + "-" + string + ")";
  979.         Dispatch->Push(*temp);
  980.         delete temp;
  981.         return Dispatch->ReturnMat();
  982. }
  983. //------------- unary minus
  984. VMatrix& operator -(VMatrix &ROp)
  985. {
  986.     ROp.Garbage("operator -");
  987.     VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  988.     for (int i = 1; i <= ROp.r;++ i) {
  989.         for (int j = 1; j <= ROp.c;++ j) {
  990.             temp->M(i, j) = -ROp.m(i, j);
  991.         }
  992.     }
  993.     temp->name = temp-> name + "-" + ROp.name + ")";
  994.     Dispatch->Push(*temp);
  995.     delete temp;
  996.     return Dispatch->ReturnMat();
  997. }
  998.  
  999. //-------------- multiplication
  1000.  
  1001. VMatrix& operator* (VMatrix &LOp, VMatrix &ROp)
  1002. {
  1003.     long double sum = 0;
  1004.     LOp.Garbage("operator M*M");
  1005.     ROp.Garbage("operator M*M");
  1006.     if (LOp.c == ROp.r) {
  1007.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  1008.         for (int i = 1; i <= LOp.r; i++) {
  1009.             for (int j = 1; j <= ROp.c; j++) {
  1010.                 sum = 0.0;
  1011.                 for (int k = 1; k <= LOp.c; k++)
  1012.                     sum +=(long double)
  1013.                         ((long double) LOp.m(i, k)) *((long double) ROp.m(k,j));
  1014.                 temp->M(i, j) = (double) sum;
  1015.             }
  1016.         }
  1017.         temp->name = temp-> name + LOp.name + "*" + ROp.name + ")";
  1018.         Dispatch->Push(*temp);
  1019.         delete temp;
  1020.     }
  1021.     else Dispatch->Nrerror(
  1022.             "Mismatched VMatrix sizes in multiplication function\n");
  1023.     return Dispatch->ReturnMat();
  1024. }
  1025. VMatrix& operator*(double scalar, VMatrix &ROp)
  1026. {
  1027.     ROp.Garbage("operator s*M");
  1028.     strtype string;
  1029.     char buffer[MAXCHARS] = { '\0' };
  1030.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  1031.         for (int i = 1; i <= ROp.r;++ i) {
  1032.             for (int j = 1; j <= ROp.c;++ j) {
  1033.                 temp->M(i, j) = scalar * ROp.m(i, j);
  1034.             }
  1035.         }
  1036.         gcvt(scalar, NDECS, buffer);
  1037.         string = buffer;
  1038.         temp->name = temp-> name + string + "*" + ROp.name + ")";
  1039.         Dispatch->Push(*temp);
  1040.         delete temp;
  1041.         return Dispatch->ReturnMat();
  1042. }
  1043. VMatrix& operator* (VMatrix &ROp, double scalar)
  1044. {
  1045.     ROp.Garbage("operator M*s");
  1046.     strtype string;
  1047.     char buffer[MAXCHARS] = { '\0' };
  1048.  
  1049.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  1050.         for (int i = 1; i <= ROp.r;++ i) {
  1051.             for (int j = 1; j <= ROp.c;++ j) {
  1052.                 temp->M(i, j) = scalar * ROp.m(i, j);
  1053.             }
  1054.         }
  1055.         gcvt(scalar, NDECS, buffer);
  1056.         string = buffer;
  1057.         temp->name = temp-> name + ROp.name + "*" + string + ")";
  1058.         Dispatch->Push(*temp);
  1059.         delete temp;
  1060.         return Dispatch->ReturnMat();
  1061. }
  1062.  
  1063. //-------- elementwise multiplication
  1064.  
  1065. VMatrix& operator %(VMatrix &LOp, VMatrix &ROp)
  1066. {
  1067.     LOp.Garbage("operator M%M");
  1068.     ROp.Garbage("operator M%M");
  1069.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  1070.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  1071.         for (int i = 1; i <= LOp.r;++ i) {
  1072.             for (int j = 1; j <= ROp.c;++ j)
  1073.                 temp->M(i, j) = LOp.m(i, j) *ROp.m(i, j);
  1074.         }
  1075.         temp->name = temp-> name + LOp.name + "%" + ROp.name + ")";
  1076.         Dispatch->Push(*temp);
  1077.         delete temp;
  1078.     }
  1079.     else Dispatch->Nrerror(
  1080.             "Mismatched Matrix sizes in elementwise multiplication\n");
  1081.     return Dispatch->ReturnMat();
  1082.  
  1083. }
  1084.  
  1085. //------------- division
  1086.  
  1087. VMatrix& operator /(VMatrix &LOp, VMatrix &ROp)
  1088. {
  1089.     LOp.Garbage("operator M/M");
  1090.     ROp.Garbage("operator M/M");
  1091.  
  1092.     if (LOp.r == ROp.r && LOp.c == ROp.c) {
  1093.         VMatrix *temp = new VMatrix("(", LOp.r, ROp.c);
  1094.         for (int i = 1; i <= LOp.r;++ i) {
  1095.             for (int j = 1; j <= ROp.c;++ j) {
  1096.                 double d = ROp.m(i, j);
  1097.                 double L = LOp.m(i, j);
  1098.                 //(fabs(d) > ZERO || fabs((d - L) / d) < 1.0E-5 )
  1099.                 temp->M(i, j) = ( 0.0 != d ) ? L / d
  1100.                               : ROp.Nrerror(" zero divisor in elementwise division");
  1101.             }
  1102.         }
  1103.         temp->name = temp-> name + LOp.name + "/" + ROp.name + ")";
  1104.         Dispatch->Push(*temp);
  1105.         delete temp;
  1106.     }
  1107.     else Dispatch->Nrerror(
  1108.             "Mismatched Matrix sizes in elementwise division\n");
  1109.     return Dispatch->ReturnMat();
  1110.  
  1111. }
  1112.  
  1113. VMatrix& operator /(VMatrix &ROp, double scalar)
  1114. {
  1115.     ROp.Garbage("operator M/s");
  1116.     strtype string;
  1117.     char buffer[MAXCHARS] = { '\0' };
  1118.  
  1119.         if (fabs(scalar) < ZERO)
  1120.             ROp.Nrerror(" zero divisor in elementwise division");
  1121.  
  1122.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  1123.         for (int i = 1; i <= ROp.r;++ i) {
  1124.             for (int j = 1; j <= ROp.c;++ j) {
  1125.                 temp->M(i, j) = ROp.m(i, j) / scalar;
  1126.             }
  1127.         }
  1128.         gcvt(scalar, NDECS, buffer);
  1129.         string = buffer;
  1130.         temp->name = temp-> name + ROp.name + "/" + string + ")";
  1131.         Dispatch->Push(*temp);
  1132.         delete temp;
  1133.         return Dispatch->ReturnMat();
  1134. }
  1135.  
  1136. VMatrix& operator /(double scalar, VMatrix &ROp)
  1137. {
  1138.     ROp.Garbage("operator s/M");
  1139.     strtype string;
  1140.     char buffer[MAXCHARS] = { '\0' };
  1141.         VMatrix *temp = new VMatrix("(", ROp.r, ROp.c);
  1142.         for (int i = 1; i <= ROp.r;++ i) {
  1143.             for (int j = 1; j <= ROp.c;++ j) {
  1144.                 temp->M(i, j) = (ROp.m(i, j) != 0.0) ? scalar / ROp.m(i, j)
  1145.                 : ROp.Nrerror("zero divisor in scalar divison");
  1146.             }
  1147.         }
  1148.         gcvt(scalar, NDECS, buffer);
  1149.         string = buffer;
  1150.         temp->name = temp-> name + string + "/" + ROp.name + ")";
  1151.         Dispatch->Push(*temp);
  1152.         delete temp;
  1153.         return Dispatch->ReturnMat();
  1154. }
  1155.  
  1156. //////////////////////end of friend functions
  1157.  
  1158. ////////////////////// matrix functions
  1159.  
  1160. ////////////// utility functions
  1161.  
  1162. void VMatrix::Writeb(char *fid, VMatrix &mat)
  1163.     /* write a matrix in an binary file */
  1164.     {
  1165.         int i;
  1166.         FILE *file;
  1167.         double *x;
  1168.         char name[MAXCHARS] = { '\0' };
  1169.  
  1170.         file = fopen(fid, "wb");
  1171.         if (!file) Dispatch->Nrerror("WRITEB: unable to open file");
  1172.  
  1173.         mat.Garbage("Writeb");
  1174.  
  1175.         strncpy(name, mat.StringAddress(), 79);
  1176.         i = strlen(name);
  1177.         fwrite(&i, sizeof (int), 1, file);
  1178.         fputs(name, file);
  1179.         i = mat.r;
  1180.         fwrite(&i, sizeof (int), 1, file);
  1181.         i = mat.c;
  1182.         fwrite(&i, sizeof (int), 1, file);
  1183.  
  1184.         unsigned blocks = mat.hdr-> nblocks;
  1185.         char *locbuf = mat.hdr-> next.v->buf;
  1186.         unsigned block_num;
  1187.  
  1188.         for (i = 0; i < blocks; i++) {
  1189.             if (!load_block(mat.hdr, (unsigned) i))
  1190.                 Nrerror("WRITEB: could not load a virtual block");
  1191.             fwrite(locbuf, sizeof (char), BLK_SIZE, file);
  1192.         }
  1193.  
  1194.         fclose(file);
  1195.         mat.M(1, 1);     // reset matrix to base pointer
  1196. }     /* writeb */
  1197.